#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "galois.h"

static int prim_poly[33] = 
{ 0, 
/*  1 */     1, 
/*  2 */    07,
/*  3 */    013,
/*  4 */    023,
/*  5 */    045,
/*  6 */    0103,
/*  7 */    0211,
/*  8 */    0435,
/*  9 */    01021,
/* 10 */    02011,
/* 11 */    04005,
/* 12 */    010123,
/* 13 */    020033,
/* 14 */    042103,
/* 15 */    0100003,
/* 16 */    0210013,
/* 17 */    0400011,
/* 18 */    01000201,
/* 19 */    02000047,
/* 20 */    04000011,
/* 21 */    010000005,
/* 22 */    020000003,
/* 23 */    040000041,
/* 24 */    0100000207,
/* 25 */    0200000011,
/* 26 */    0400000107,
/* 27 */    01000000047,
/* 28 */    02000000011,
/* 29 */    04000000005,
/* 30 */    010040000007,
/* 31 */    020000000011, 
/* 32 */    00020000007 };  /* Really 40020000007, but we're omitting the high order bit */

main (int argc, char **argv)
{
	int remainder[2][5], auxiliary[2][5];

	int w;
	int w32;
	int a;
	int high_bit;
	int high_bit_r;
	int r1, q1, shift_count, a1;

	if (argc != 2) { fprintf(stderr, "usage: prog w\n"); exit(1); }

	w = atoi(argv[1]);
	if (w > 32) exit(1);
	if (w == 32) w32 = 1; else w32 = 0;

	for(a = 2; a < (1 << w); a++) {
		remainder[0] = prim_poly[w];
		remainder[1] = a;

		auxiliary[0] = 0;
		auxiliary[1] = 1;

		high_bit = 4;
		high_bit_r = 4;
		while (remainder[1] > 1) {
			q1 = 0;
			high_bit = high_bit_r;

			for (high_bit_r = 3; remainder[1][high_bit_r] == 0; high_bit_r--);

			/* This is the polynomial division part */
			while(high_bit >= high_bit_r) {
				r1 = remainder[1];
				shift_count = 0;

				/* If w == 32, the width is bigger than an int, so deal with that: */
				if (w32) {
					while((r1 & (1 << (high_bit-1))) == 0) { r1 <<= 1; shift_count++; }
					r1 <<= 1; shift_count++;
					w32 = 0;
				} else while((r1 & (1 << high_bit)) == 0) { r1 <<= 1; shift_count++; }

				/* set the quotient */
				q1 ^= (1 << shift_count);

				/* next remainder in the division process */
				remainder[0] ^= r1;
				while((remainder[0] & (1 << high_bit)) == 0) high_bit--;
			}

			r1 = remainder[0];
			remainder[0] = remainder[1];
			remainder[1] = r1;
			a1 = galois_single_multiply(q1, auxiliary[1], w);
			a1 ^= auxiliary[0];
			auxiliary[0] = auxiliary[1];
			auxiliary[1] = a1;
		}
		printf("%3d with inverse %3d\n", a, auxiliary[1]);
	}
}
